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Abstract. In this paper we present a numerical algorithm for solving 
the direct scattering problems by the Modified Rayleigh Conjecture Method 
(MRC) introduced in [1]. Some numerical examples are given. They show 
that the method is numerically efficient. 
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I. Introduction 

The classical Rayleigh Conjecture is discussed in [4] and [5], where it is 
shown that, in general, this conjecture is incorrect: there are obstacles (for 
example, sufficiently elongated ellipsoids) for which the series, representing 
the scattered field outside a ball containing the obstacle, does not converge 
up to the boundary of this obstacle. 



The Modified Rayleigh Conjecture (MRC) has been formulated and 
proved in [1] (see Theorem 1 below) . A numerical method for solving obsta- 
cle scattering problems, based on MRC, was proposed in [1]. This method 
was implemented in [2] for two-dimensional obstacle scattering problems. 
The numerical results in [2] were quite encouraging: they show that the 
method is efficient, economical, and is quite competitive compared with 
the usual boundary integral equations method (BIEM). A recent paper [3] 
contains a numerical implementation of MRC in some three-dimensional ob- 
stacle scattering problems. Its results reconfirm the practical efficiency of 
the MRC method. 

In this paper a numerical implementation of the Modified Rayleigh Con- 
jecture (MRC) method for solving obstacle scattering problem in three - 
dimensional case is presented. Our aim is to consider more general than in 
[3] three-dimensional obstacles: non-convex, non-starshaped, non-smooth, 
and to study the performance of the MRC in these cases. The minimiza- 
tion problem (5) (see below), which is at the heart of the MRC method, is 
treated numerically in a new way, different from the one used in [2] and [3] . 
Our results present further numerical evidence of the practical efficiency of 
the MRC method for solving obstacle scattering problems. 

The obstacle scattering problems (l)-(3), we are interested in, consists 
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of solving the equation 

{^^ + k'^)u = in D' = R^\D, (1) 

where D C is a bounded domain, satisfies the Dirichlet boundary condi- 
tion 

u\s = 0, (2) 

where S is the boundary of D, which is assumed Lipschitz in this paper, 
and the radiation condition at infinity: 

u = uq + V = uo + A{a' ,a) \- o{-) r — oo, (3) 

r := \x\, a' = x/r, uq := e*'^"'^, 

where v is the scattered field, a G is given, is the unit sphere in R^, 
k = const > is fixed, k is the wave number. The coefficient A{a', a) is 
called the scattering amplitude. 
Denote 

Ma):= [ Aia',a)Yii^da', (4) 
where Yi{a) are the orthonormal spherical harmonics: 

Yi = Yi^, -l<m<l, / = 0,1,2,... 

V ZTT 



Pl^{x) are the associated Legendre functions of the first kind, 
Prix):=il-x'r^'^^^, m>0, 

and 

For m < 

ei^(a;) = (-l)-e^,_„(a;). 

Let hi{r) be the spherical Hankel functions of the first kind, normahzed so 
that hi{kr) ~ e^'^^/r as r ^ +00. Let Br := {x : \x\ < R} D D, and the 
origin is inside D. 

Then in the region r > R, the solution to the acoustic wave problem 
(l)-(3) is of the form: 

00 

u{x,a) = e^'^"'^ + '^Ai{a)'il}i{x), \x\ > R, 
1=0 

ipl := Yi{a')hi{kr), r > R, a' = x/r, 

where 

00 00 z 

E-E E • 

1=0 1=0 m=-l 

Fix e > 0, an arbitrary small number. The following Lemmas and Theoreml 
are proved in [1]. 



Lemma 1. There exist L = L{e) and numbers ci = q(e) such that 

L 

Iko(s) + X]cKe)'0Ks)IU2(5) <e. (5) 

1=0 

Lemma 2. If (5) holds, then \\ve{x) - v{x)\\ = 0{e), \/x e D', e ^ 0. 

where 

L 

:=^Q(e)V'/(a;), xeD', (6) 

1=0 

and 

IMI := IMIh^'^^P') + IMlL2(D';(l+|a;|)-^)' 7 > 0, m > 0, (7) 

m is arbitrary, and H"^ is the Sobolev space. 
Lemma 3. c;(e) Ai{a),\/l,e 0. 

Theorem 1 (Modified Rayleigh Conjecture). Let D e be a 

bounded obstacle with Lipschitz boundary S. For any e > there exists L = 
L{e) and q(e) = Q^(e), Q < I < L, —I < m < I, such that inequality (5) 
holds. If (5) holds then function (6) satisfies the estimate \ \v{x) — v^{x)\ \ = 
0{e), where the norm is defined in (7). Thus, v^{x) is an approximation of 
the scttered field everywhere in D' . 

In order to obtain an accurate solution, usually one has to take L large. 
But as L grows the condition number of the matrix (^i, tpi')L^(s) is increas- 
ing very fast. So we choose some interior points xj e D, j = 1,2, ...,J, and 
use the following version of Theorem 1([2]): 
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Theorem 2. Suppose xj E D, j = 1,2, J, then Ve > 0, 3L = L(e) 
and cij{e),l = 0, L, j = 0, J(e), such that 
(i) 

J L 

IKo(s) + X!Z!^o(«)^^(*~^i)lli2(5) <e- (5') 
|K(x)-i;(x)|| = 0(e), 

where 

J L 
j=Q 1=0 

and the \ \.\\ is defined in Lemma 2. 

Remark. Theorem 1 is the basis for MRC algorithm for computation 
of the field scattered by an obstacle: one takes an e > and an integer 
L > 0, minimizes the left-hand side of (5) with respect to q, and if the 
minimum is < e then the function (6) is the approximate solution of the 
obstacle scattering problem with the accuracy 0(e). If the above minimum 
is greater than e, then one increases L until the minimum is less than e. 
This is possible by Lemma 1. In computational practice, one may increase 
also the number J of points Xj inside D, as explained in Theorem 2. The 
increase of J allows one to reach the desired value of the above minimum 
keeping L relatively small. This gives computational advantage in many 
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cases. 



In section 2, an algorithm is presented for solving the problem (l)-(3). 
This algorithm is based on the MRC. Compared with the previous work 
in the case of two- and three-dimensional MRC ([2], [3]), we consider more 
general surfaces, in particular non-starshaped and piecewise-smooth bound- 
aries. The numerical results are given in section 3. A discussion of the 
numerical results is given in section 4. 



II. The MRC algorithm for Solving Obstacle Scat- 
tering Problems 



1. Smooth starshaped boundary: 



Assume the surface S is given by the equation 



r = r{9,ip) 



< (ys < 27r, < 6" < TT. 



Define 



L 




(5") 



Let 



hi = 27r/ni, 



h2 = 7r/n2 



= (pQ<(pi < .... < (fini = 27r, (fii^ = iihi, ii = 1, ...,ni, 



= 0Q < 9l < .... < 6ni = TT, = i2h2, 12 = 1, 
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where ni and n2 are the number of steps. By Simpson's formula([8]), we 
obtain an approximation of F{co, ci, c^): 



ni n2 



F(co,ci, ...,cl) = XI XI "nj2l'"0ni2 + X '^^^fe»2l^^nj2^i^2 (5"') 

ii=0i2=0 1=0 



where 



and 



4, ii and 12 even 
Oii.io = < 8, h — i2 odd 
16, zi and 12 odd 



where 



w{9, (p) = {r'^rl + r'^rlsin^e + r'^sin^ef/'^ . 



We can find c* = {cq,c\, ...,c\) such that 



(8) 



F{c*) = minF{cQ, ci, cl). 



(9) 



We first write 



F(c) = \\Ac-B\\^, 



(10) 



where 



A = {Ai^i)MxLi, Ai^i = ipii^i2{ai^i^Wi^i^hih2)2 , i = 1112, 



B = {Bi)Mxii Bi = uoi^i^{ai^i^Wi^i^hih2)'^ , 

in which M = n\n2, Li = {L + 1)(2L + 1) since q = q^, < Z < L, — Z < 
m < I. 

Then Householder reflections are used to compute an orthogonal-triangular 
factorization: A * P = Q * R where P is a permutation([8], p. 171), Q 
is an orthogonal matrix, and R is an upper triangular matrix. Let r = 
rank{A). This algorithm requires 4MLir - 2r'^{M + Li) + 4r2/3 flops([9], 
pp. 248-250). The least squares solution c is computed by the formula 
c = P * {R~^ * {Q' * {A^B))). This minimization procedure is based on 
the matlab code([10]). 

In [2] and [3] singular value decomposition was used for minimization 
of (5"). Here we use the matlab minimization code which is based on a 
factorization of the matrix A. This has the following advantages from the 
point of view of numerical analysis. We can choose an integer ri: 

< ri < r 

such that the first ri rows and columns of R form a well-conditioned matrix 
when A is not of full rank, or the rank of A is in doubt([10]). See Golub 

and Van Loan [9] for a further discussion of numerical rank determination. 
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If we choose xj G D, j = 1, J, we obtain 

Fj{c) = Fj{coi, ...,coj,cn, ■■■,cij, ...,cli, ■■■,clj) 

n\ n2 J L 
ii=0i2=0j=l 1=0 

The algorithm for finding the minimum of Fj(c) will be same. 

2. Piecewise-smooth boundary: 
Suppose 

N 

S=\JSn. 

n=l 

Then 

AT L 

F{co,ci,...,cl) = J2\\^o + J2^^i^i\\L^Sn) 

n=l 1=0 

y{x,y,z)ESn, = + + , cos9 = z/r, tan(p = y/x. (11) 

3. Non-starshaped case: 

Suppose 5 is a finite union of the surfaces, each of which is starshaped 
with respect to a point r^, 

JV 

^ = U 

n=l 

and the the surfaces Sn are given by the equations in local spherical coordi- 
nates: 



H) ¥'n) cos c/?7i sin ^jj, r7i(^fj, sin (^jj sin ^jj, rn(On->'^n) ^^^(^n)t 

n = 1,...N, 
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where are constant vectors. 
Then 

N L 

F{co,ci,...,cl) = \ \uo + J2'^i^i\\Ihs^)- 

n=l 1=0 

The weight functions Wn{0, ffi) are the same as in (8) since r° are constant 
vectors. 

III. Numerical Results 

In this section, we give four examples to show the convergence rate of 
the algorithm and how the error depends on the shape of S. 

Example 1. The boundary S is the sphere of radius 1 centered at the 
origin. 

In this example, the exact coefficients are: 

Clm = , n \ '^ImK.Oi) 

hi(k) 

Let A; = 1, a = (1, 0, 0). We choose ni = 20, n2 = 10. 



L 


1.0000 


2.0000 


3.0000 


4.0000 


5.0000 


6.0000 


7.0000 


F{c*) 


6.3219 1.6547 


0.2785 


0.0368 


0.0034 


0.0003 


0.0000 


0.0000 



err(c) 0.0303 0.0172 0.0020 0.0004 0.0000 0.0000 0.0000 0.0000 
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where 

L ^ 

err(c) = (X^lcf -Qp)^. 

1=0 

When ni = 40, n2 = 20, 



L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 



F{c*) 6.3544 1.6562 0.2820 0.0358 0.0036 0.0003 0.0000 0.0000 



err(c) 0.0147 0.0076 0.0011 0.0001 0.0000 0.0000 0.0000 0.0000 



Next, we fix ni = 20, n2 = 10 and test the results for different k and a. 
When A; = 2, a = (1,0,0), 



L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 



F{c*) 10.4506 5.5783 1.9291 0.5217 0.0970 0.0156 0.0020 0.0003 



err(c) 0.0404 0.0205 0.0048 0.0020 0.0005 0.0000 0.0000 0.0000 
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When A; = 1, a = (0,1,0), 



L 





1.0000 2.0000 


3.0000 


4.0000 


5.0000 


6.0000 


7.0000 


F(c*) 


6.3801 


1.6628 0.2821 


0.0371 


0.0044 


0.0003 


0.0000 


0.0000 


err(c) 


0.0014 


0.0106 0.0005 


0.0004 


0.0000 


0.0000 


0.0000 


0.0000 


When A; = 1, 


« = (0,0,1), 












L 





1.0000 2.0000 


3.0000 


4.0000 


5.0000 


6.0000 


7.0000 


F{c*) 


6.4156 


1.6909 0.2955 


0.0418 


0.0025 


0.0002 


0.0000 


0.0000 


err(c) 


0.0093 


0.0109 0.0049 


0.0007 


0.0001 


0.0000 


0.0000 


0.0000 


When A; = 1, 


a = (1/^2, 1/V2, 0), 










L 





1.0000 2.0000 


3.0000 


4.0000 


5.0000 


6.0000 


7.0000 
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F{c*) 6.3500 1.6711 0.2810 0.0371 0.0040 0.0003 0.0000 0.0000 



err(c) 0.0218 0.0057 0.0019 0.0004 0.0001 0.0000 0.0000 0.0000 



When A; = 1, 


a = (1/^3,1/^3,1/^3), 








L 


1.0000 2.0000 3.0000 4.0000 


5.0000 


6.0000 


7.0000 


F{c*) 6.3739 


1.6542 0.2850 0.0368 0.0040 


0.0003 


0.0000 


0.0000 


err(c) 0.0170 


0.0054 0.0021 0.0003 0.0001 


0.0000 


0.0000 


0.0000 



Example 2. The boundary S is the surface of the cube [—1, 1]^. Here 

6 

S=\JSn. 
n=l 

and 

6 L 

F{co,ci,...,cl) = W^O + ^Ci^lWh^Sn) 
n=l 1=0 

6 ni H'z L 
n=lii=0i2=0 1=0 

where 

Ai = 2/ni, A2 = 2/n2. 
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The origin is chosen at the cenetr of symmetry of the cube. The surface area 
element is calculated in the Cartesian coordinates, so the weight w = 1. 
Let be the surface 

z = l, -l<a;<l, -l<y<l 

and 

= — 1 + iiAi, < ii < rii 
yi2 = -1 + i2^2, < i2 < n2 

Then 

V'Hiia = Yl{Gh,Vi2)h{kr{ei^,ipi^)), 

and and (pi^ can be computed by formula (11). For other surfaces Sj the 
algorithm is similar. 

The values of mini<'(c) = F{c*) and the values minFj(c) = Fj[c*) with 

Xj: 

{x^ : j = 0, 6} = {(0, 0, 0), (0.2, 0, 0), (-0.2, 0, 0), 

(0, 0.2, 0), (0, -0.2, 0), (0, 0, 0.2), (0, 0, -0.2)} 

are given below. 

We choose ni = 10, n2 = 10 
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L 







1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 



F{c*) 10.6301 3.6277 2.6760 2.2309 1.9832 1.5737 1.5034 1.2948 1.1753 



Fj{c*) 2.6297 1.0970 0.5487 0.1572 0.0667 0.0320 0.0168 0.0078 0.0035 



When m = 20, n2 = 20, 



L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 



F{c*) 10.7923 3.7144 2.7778 2.3393 2.0873 1.6671 1.5938 1.4277 1.3368 



Fj{c*) 2.7248 1.1433 0.5757 0.1686 0.0694 0.0652 0.0236 0.0143 0.0090 



Example 3. The boundary S is the surface of the eUipsoid x"^ + y"^ + 
z^/b"^ = 1, the values of minF(c) = -F(c*), 6 = 2, 3,4, 5 with m = 20, n2 = 
10 are: 



L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 
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b=2 8.8836 5.4955 3.0421 2.8434 1.3622 1.2093 0.8753 0.8132 



b=3 14.1617 12.0477 7.2296 7.0999 3.8077 3.6829 3.1324 3.0496 



b=4 19.5326 17.9346 9.9927 9.8720 5.3333 5.2008 4.6793 4.5738 



b=5 22.9765 21.5653 11.4850 11.3587 6.1637 6.0096 5.5202 5.3933 



The values of min Fj(c) = Fj{c*), b = 2,3,4,5 with xj: 

{xj : i = 0, 6} = {(0, 0, 0), (0.5, 0, 0), (-0.5, 0, 0), 
(0, 0.5, 0), (0, -0.5, 0), (0, 0, 0.5), (0, 0, -0.5)} 

are: 



L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 



b=2 2.4856 0.7090 0.2530 0.0062 0.0000 0.0000 0.0000 0.0000 
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b=3 4.6639 1.3619 0.6618 0.0074 0.0000 0.0000 0.0000 0.0000 



b=4 5.5183 1.8624 0.7844 0.0060 0.0000 0.0000 0.0000 0.0000 



b=5 11.0579 8.7027 6.4831 0.8357 0.0017 0.0000 0.0000 0.0000 



Example 4. The obstacle is a dumbbell. Its boundary S is not smooth, 
non-starshaped and not convex: 

S = Si\JS2\JS3 

Si : r — (0, 0, 1) = (1.5 cos ipsin 9, 1.5 sin sin ^, 1.5 cos 6) 
82'- f — (0, 0, —1) = (1.5cos </3sin^, l.^snnpsmO, 1.5 cos 0) 

: rsinO = 1 

{xj : i = 0, 10} = {(0, 0, 0), (0, 0, 0.1), (0, 0, -0.1), (0, 0, 0.2), (0, 0, -0.2), 

(0, 0, 0.3), (0, 0, -0.3), (0, 0, 0.4), (0, 0, -0.4), (0, 0, 0.5), (0, 0, -0.5)}; 
We choose ni = 20, n2 = 10 for every = 1, 2, 3). 
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L 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 



F{c*) 


25.8840 


20.8059 


16.4968 


15.6622 12.9241 12.1915 11.0187 9.5263 


Fj{c*) 


20.3118 


8.0238 


5.1062 


2.5908 0.8304 0.4067 0.0453 0.0084 



IV. Conclusion 

Prom the numerical results one can see that the accuracy of the numerical 
solution depends on the smoothness and elongation of the object. 

In Example 1 the surface 5 is a unit sphere and the numerical solution 
is very accurate. In Example 3 the results for different elongated ellipsoids 
show that if the elongation (eccentricity) grows, then the accuracy decreases. 
In Example 2 the surface is not smooth and the result is less accurate than 
in Example 3. In Example 4 the surface in nonconvex and not smooth, but 
the accuracy is of the same order as in Example 2. 

When b is large or S is not smooth, the numerical results in Example 2 
and Example 3 show that if one adds more points Xj then the accuracy of 
the solution increases. 

In Example 1 and Example 2, as one increased ni and n2, the minimum 
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F{c*) has also increased because the condition number of the matrix A in 
(10) grew as ni and n2 increased. 

Using the results of Example 1 one can check the accuracy in finding c/ 
by the value of the minimum 

F{c*) < e. 
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